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ABSTRACT 

Context. Disk accretion onto weakly magnetized stars leads to the formation of a boundary layer (BL) where the gas loses its excess 
kinetic energy and settles onto the star. There are still many open questions concerning the BL, for instance the transport of angular 
momentum (AM) or the vertical structure. 

Aims. It is the aim of this work to investigate the AM transport in the BL where the magneto-rotational instability (MRI) is not 
operating owing to the increasing angular velocity Dfr) with radius. We will therefore search for an appropriate mechanism and 
examine its efficiency and implications. 

Methods. We perform 2D numerical hydrodynamical simulations in a cylindrical coordinate system (r, (p) for a thin, vertically inte¬ 
grated accretion disk around a young star. We employ a realistic equation of state and include both cooling from the disk surfaces and 
radiation transport in radial and azimuthal direction. The viscosity in the disk is treated by the a-model; in the BL there is no viscosity 
term included. 

Results. We find that our setup is unstable to the sonic instability which sets in shortly after the simulations have been started. Acoustic 
waves are generated and traverse the domain, developing weak shocks in the vicinity of the BL. Furthermore, the system undergoes 
recurrent outbursts where the activity in the disk increases strongly. The instability and the waves do not die out for over 2000 orbits. 
Conclusions. There is indeed a purely hydrodynamical mechanism that enables AM transport in the BL. It is efficient and wave 
mediated; however, this renders it a non-local transport method, which means that models of a effective local viscosity like the 
a-viscosity are probably not applicable in the BL. A variety of further implications of the non-local AM transport are discussed. 
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1. Introduction 

One ubiquitous phenomenon in astrophysics is the accretion of 
matter on a central object via an accretion disk. This process can 
be observed for a variety of central objects, such as young stars, 
compact objects like white dwarfs or neutron stars, and active 
galactic nuclei. The physics of the accretion disk itself is reason¬ 
ably well understood. However, the tiny region where the accre¬ 
tion disk connects to the star is still one of the major problems. In 
this region, called the boundary layer (BL), a smooth connection 
is established between the disk, which rotates nearly at Keple- 
rian frequency, and the star, which in general rotates much more 
slowly. It is of great dynamical and thermal importance since the 
gas undergoes a supersonic velocity drop in a very confined re¬ 
gion of a few percent of the stellar radius and loses up to one 


half of the total accretion energy during that process (Kluzniak 
T987l l. 

The BL has been studied extensively for over 40 years, both 
in an analytic al (e.g. |Bertout & Regev|1992[ ) and a numerical 
approach (e.g. Kley & Hensler|1987| l. Simulations involving the 
BL around a young star have been performed by |Kley & Lin| 
(|1996|l and|Kley & Lin|(|1999[), including full radiation transport. 


Most of the work that has been done on the BL (see e.g. Hert 


|felder et al.|2013| for a brief review of the history) assumes that 
the gas first slows down in the equatorial plane of the disk and 
then spreads around the star. This is the classical picture of the 
BL. However, there is an alternative theory called the spread¬ 


ing layer ( jlnogamov & Sunyaev| 1999] 20101, which was ini¬ 
tially developed for neutron stars. Within this concept the gas 
first spreads around the star because of the ram pressure in the 
BL and then slows down on the whole surface of the star. 

The majority of the simulations of the BL start from the 
premise that the observed angular momentum (hereafter AM) 
transport in the BL is driven by local turbulent stresses and 
the authors consequently adopt a a-viscosity model ( |Shakura'&] 
Sunyaev|1973|l in their simulations. Sometimes the classical a- 


model is modified slightly in order to prevent supersonic infall 
velocities ( |Kley & Papaloizou| 1997) 1 or, for instance, to take into 
account the radial pressure scale height in the BL (|Popham'& 


|Narayan| 19951. The utilization of a local viscosity model was 


later justified for the disk by the discovery of the magnetorota- 
tional instability (MRI) which creates turbulence that acts like a 
genuine viscosity on macroscopic scales (|Velikhov|1959HChan-| 
drasekhar|[T960HBalbus & Hawley||1991| |1998| |Balbus||2003[ 


Balbus & Lesaffre 2008|l. However, as pointed out by 


(1995 I and Abramowicz et al.|(|1996|l and recently shown by|Pes 


Godon 


|sah & Chan|(|2012|l, if the angular velocity increases with radius, 
dG/dr > 0, the MRI is effectively damped out and the asso¬ 
ciated AM transport oscillates about zero. Since this situation 
clearly applies for the BL, we do not expect to obtain sufficient 
AM transport through the MRI. There have been various alter¬ 
native transport mechanisms proposed, among them the Kelvin- 
Helmholtz instability (|Kippenhahn & Thomas||1978|l, the baro- 


clinic instability (Lujimoto 
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1993|), and the Tayler-Spruit dynamo 


( |Tayler| 19^3)[Spruit 2002[ Piro & Bildsten|26(j4) l, none of which 
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have yet been proven to efficiently transport mass and AM in the 
BL. 

In this work, we will focus on the AM transport in the BL 
through non-axisymmetric instabilities and therefore carry out 
2D simulations in the midplane of the star-disk system using 
cylindrical coordinates. Recently, a promising candidate for the 
transport has been proposed and investigated in a series of pa¬ 
pers. According to this theory, the steep velocity drop in the 
BL is prone to the sonic instability (Glatzel |1988[ [Belyaev ^ 


|Rafikov||2012| ), which is an instability of a supersonic shear 
layer, much like the Papaloizou-Pringle instability (jPapaloiz^ 


& Pringle||1984| Narayan et al.||1987) l. Acoustic waves are ex¬ 
cited in the BL as a consequence of the sonic instability and 
AM can be transported by these modes in an efficient way. This 
has been demonstrated for 2D flows (Belyaev et al. 2012| l, 3D 
flows in cylindrical coordinates ( [Belyaev et al.|2013a I, and even 
for 3D magnetohydrodynamical flows (Belyaev et al. 2013b[ l. 
The wave mediated AM transport implies that this process is 
intrinsically non-local since the waves can potentially travel a 
long way before they dissipate and release the AM to the fluid. 
Therefore, it is problematic to describe the AM transport in the 
BL by means of a local viscosity like the a-model. Although the 
above-mentioned simulations are already quite sophisticated, the 
authors make some simplifications that constrain the validity of 
their models. This is the point where we step in with this paper 
and relax three of the simplifications made. First, we make use 
of a realistic equation of state instead of an isotfiermal one and 
propagate tfie temperature of tfie system, as well. We use full ra¬ 
diative diffusion in tfie disk plane and an approximation for the 
vertical flux, and hence employ a quasi-3D radiation transport. 
Second, we use realistic, state-of-the-art initial models like those 
already published in [Hertfelder et ak ( 2013| l as a starting point 
for the simulations presented here. Third, mass constantly enters 
the simulation domain from the outer boundary and thus there is 
a net mass flow through the disk, as is the case for a real accre¬ 
tion disk. Because of the high computational costs the treatment 
of the radiation brings with it, we perform 2D simulations and 
do not yet include magnetic fields. Thus, the next steps will be to 
extend the simulations presented here to 3D and also to include 
magnetic fields. 

The paper is organized as follows. In Sect. 2, the basic phys¬ 
ical model is described. We give an overview of the equations 
used for our simulations and present the assumptions we have 
made. Furthermore, we briefly review how the AM transport due 
to stresses in the fluid can be measured and quantized. Section 3 
is devoted to the numerical method that we employ in our simu¬ 
lations. We discuss the numerical code, the initial and boundary 
conditions, as well as the model parameters. In Sect. 4 we start 
the analysis of the results. We first discuss tfie sonic instabil¬ 
ity and the acoustic modes that are excited in the BL. Then, in 
Sect. 5, we investigate how these instabilities affect the physi¬ 
cal picture of the BL and look deeply into the AM transport and 
long-term behavior of our models. We conclude with Sect. 5. 


2. Physics 

In this section, we present the physical foundations for the sim¬ 
ulations we have performed and that are described later in this 
publication. In order to investigate non-axisymmetric behavior 
in the disk plane, we utilized the vertically integrated Navier- 
Stokes equations. Since the focus of our current work lies on 
instability developing in the BL, we describe the flow of a per¬ 
turbed, compressible fluid and also the derivation of the resulting 
Reynolds stresses. 
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2.1. Two-dimensional equations fora flat disk 

In order to obtain a set of 2D Navier-Stokes equations includ¬ 
ing radiation transport, we apply a cylindrical coordinate system 
with coordinates (r, (f, z) and integrate the 3D equations over the 
vertical coordinate z. The star lies in the center of the coordinate 
system and the midplane of the disk coincides with z = 0. If we 
assume a Gaussian profile for the mass density p, we can define 
the 2D surface density by 

X CXD 

pdz- "'fhtp(r,ip,0)H, (1) 

OO 

where we introduce the effective pressure scale height H, which 
is a measure of the vertical extent of the disk. Under the assump¬ 
tion of hydrostatic balance and an isothermal equation of state in 
the z-direction, the pressure scale height reads 



( 2 ) 


with Cs and Gk = s/GmTJi^ being the midplane soundspeed 
and the Keplerian angular velocity, respectively. The inverse of 
Gk(^*) is proportional to the period of a Keplerian orbit at the 
surface of the star, P* = InjTlYXRf). Derivatives with respect to 
z are neglected, as is the z-component of the velocity, u,. Clearly, 
this is not a reasonable assumption if one wants to investigate the 
overall structure of the BL. However, it is still justified since we 
restrict ourselves to tfie study of non-axisymmetric flow structure 
in tfiis work. 

The vertically integrated continuity equation reads 

(92 1 d{rl.ur) 1 

dt r dr r dtp 

We express the momenta equations in terms of the conservative 
variables and introduce the radial momentum density s = 'Lur 
and the angular momentum density h - l-ru^ = Sr^Q. These 
variables reflect the physically conserved quantities, as opposed 
to the primitive variables Ur and The conservation of momen¬ 
tum in radial direction is then given by 


(9s ^ 1 dfrsur) ^ 1 d(suf) ^u'^ Qp 
dt r dr r d<p r dr 


1 direr rr) 1 dicT rf) 

r dr r dip 


_ 1 


(4) 


and the equation in azimuthal direction reads 


dh 1 dirhuf) 1 dihuf) dp 

dt r dr r dp dp 

1 dir^a-rf) dicr^f) (90 

Here p denotes the vertically integrated pressure p - "LRcTjp, 
where Rq - Iniu with Boltzmann’s constant ks and the mass 
of hydrogen mu, T is the temperature, and p is the mean molec¬ 
ular weight; O = -GM^jr is the gravitational potential (G and 
M* are the gravitational constant and stellar mass). The compo¬ 
nents of the vertically integrated viscous stress tensor, cr,y, can be 
found in e.g. [Mass^([2002|l. The equation for the conservation 
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of energy is given by 


de 1 d(reur) 

~dt 7 dr 


+ a-y. 


dur 

"Ih 


+ cr, 


np 


1 d(eu^) 


1 drUy 

1 du^ 

r dip 

P 

r dr 

r dip 

1 duy do. 
r dip dr 

"1" O' 

do. Uy 

dip r 

- 2cr srT 

4 

d(rFy) 


dr 


dF^ 

dip 


( 6 ) 


where ctsb and Tefr are the Stefan-Boltzmann constant and the 
effective temperature, and the factor H follows from the vertical 
integration. 

The last line of Eq. (|^ incorporates the radiative flux in the 
flux-limited diffusion approximation, 


F = -K^VT. 


(7) 


Utilizing this approach, we assume that the BL is optically thick 
both in the radial and vertical directions, where the inclusion 
of a flux-limiter also allows for the treatment of optically thin 
regions. Equation Q has exactly the same mathematical form 
as molecular heat conduction in a gas. The effective radiative 
conductivity is 




4AacT^ 

krP 


( 8 ) 


where kr, a, and c are the Rosseland mean opacity, the radiation 
constant and the speed of light, and A is the flux limiter ( [Lever- 1 
more & Pomranin^ 1981 Levermor^j 1 984|l . We adopt the for- 


mulation given by jKleyj ( 1989| l for T We assume that the disk 
locally radiates like a blackbody of temperature Teff and asso¬ 
ciate the vertical flux with the blackbody luminosity divided by 
the radiating area, - crsBT^g. The factor 2 in Eq. (^ stems 
from the fact that the disk has two sides. The vertical mx pro¬ 
vides the cooling of the disk, while the other parts of the energy 
equation generate heat through viscous dissipation (second line 
in Eq. or redistribute it in the disk plane (the terms with Fy 
and F^p). The effective temperature Teff is evaluated according to 
|Hubeny| ( |1990| l, who approximates the optical depth in the ver¬ 
tical direction of the disk (see Sect. |5.3| ). By approximating the 
flux in the vertical direction as well, we employ a quasi-3D radi¬ 
ation transport. 

We also adopt a piecewise prescription for the Rosseland 
mean opacity, where in each temperature regime a-r is gi ven by 
a power-law approximation, kr = kqp^T^’, as described in Lin & 


|Papaloizou] ( [1985| ) and [Bell & Lin| ( |T994| l. Each of the regimes 
is characterized by a prevailing mechanism that dominates the 
opacity in this region. The parameters kq, a, and b for the various 
opacity regimes are listed in e.g.jMiiller & Kle^(|2012|). 


2.2. Viscosity 

In the simulations presented here, we strictly distinguish be¬ 
tween the BL and the disk with regard to viscosity. As has al¬ 
ready been mentioned in Sect. 1, the AM transport in the disk 
is mainly driven by turbulent stresses due to the MRI and can 
therefore be parametrized with the classic a-viscosity approach 
by jShakura & Sunyaev| ( |19'73| l. The effective kinematic viscosity 
is written as 

V = aCsH (9) 

in this ansatz, where a is a dimensionless parameter whose value 
is derived from MRI simulations. The cause of the AM in the BL 


and its investigation is, however, the subject of this work. Con¬ 
sequently, we do not make use of an additional viscosity term in 
the BL region. The entire AM transport in the BL is therefore 
due solely to the Reynolds stresses exerted by the gas. 

Erom a numerical point of view, the distinction between the 
BL and the disk is made by searching for the global maximum 
of the radial angular velocity profile, Q.(r,n^) = flmax- This has 
to be done for every angle ip in order to prevent a suppression 
of non-axisymmetric perturbations. We can then computation¬ 
ally realize the viscosity as pointed out above by utilizing the 
following conditions: 

- If r < Tniaxj we Set a = 0 in Eq. 0 and no viscosity is applied 
(BL region). 

- If r > Tniax, a is set to the default value, a - 0.01, and the 
viscous stress tensor does not vanish (disk region). 


2.3. Reynoids stresses 

Gas flows in general can be characterized by a laminar bulk flow 
with superimposed perturbations. One typical approach is to de¬ 
compose the flow variables / into an averaged or filtered part / 
and a fluctuating or sub filtered, unresolved part /': 

f(r,t)^f(r,t)+f(r,t). (10) 


Within the classical Reynolds approach, the averaging is per¬ 
formed over an ensemble of systems, thus / is an ensemble 
mean. In the case of a statistically steady or stationary flow field, 
or when it is homogeneous, ergodicity can be assumed. This 
means that both a sufficiently long time average over many time 
scales and a spatial average over many length scales correspond 
to the ensemble mean. Often, practical flows are inhomogeneous 
and the ensemble mean is replaced by a time average: 

/r(r,f;T)=- r f(r,t’)At’ . (11) 

T Jf-T 


If, on the other hand, there is a symmetry in at least one direction, 
as there is in our case, a spatial average is also a viable solution. 

The actual average process in the Reynolds decomposition is 
a simple arithmetic mean over space or time (see Eq.[TT]l. A more 
common approach in the case of compressible flows, however, is 
to use a density-weighted Reynolds average for the filtered part 
of the variables (jGatski & Bonnet 2013[ Peyret 1996 Balbus 


|& Papaloizou|1999| l. These density-weighted, or Favre variables 
are given by a decomposition, which reads 


f^f+r with 


p 


( 12 ) 


where the bar denotes a Reynolds average and the tilde a Eavre 
average. 

After performing a decomposition (of any kind) of the flow 
field, one can then derive the Navier-Stokes equations for the 
mean or resolved field. They are called the Reynolds-averaged 
Navier-Stokes (RANS) equations, or, in the case of a Eavre 
mean, the Eavre-averaged Navier-Stokes (EANS) equations (e.g. 
jGatski & Bonnet [201 3| |Peyret||1996| l. Those mean field equa¬ 
tions contain additional terms, however, that are dependent on 
the fluctuating part of the flow variables. In many engineering 
applications where simulations of the mean field are performed, 
those parts are unknown and have to be dealt with by using ap¬ 
propriate closure relations or models. Here, on the other hand, 
we are running direct numerical simulations and solve for the 
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actual values of the flow field variables. From these simulations, 
we can compute the additional terms, which appear in the aver¬ 
aged equations. 

From the FANS equations, we obtain the Reynolds stresses 
exerted by the fluctuating motion of the gas. For accretion disks, 
the dominant contribution is given by 


f^rtp — , (13) 

where we have already performed a vertical integration. The 
quantity R^p appears in the angular momentum equation p i, 
where it adds to cr^^. We note that for the second part in Eq. (fO i, 
it is incorrect to use either the mean values of the laminar flow 
(i.e. the initial values, when starting from a relaxed model) or 
assume Ur ~ 0. Our simulations showed that the instability- 
dominated state has other mean values than the laminar flow 
state, and that both terms in Eq. 0 are approximately of the 
same order of magnitude. Our simulations also showed that for 
the quasi-stationary state both temporal and azimuthal averaging 
procedures yield the same mean values since the evolution of the 
BE is significantly slower than an orbital period. Therefore, we 
mostly used azimuthal averaging for the mean values since it is 
computationally far less expensive. 

The vertically integrated Reynolds stress, which has the di¬ 
mensions of a 2D pressure, can be expressed by a dimensionless 
parameter aRe via 


ttRe 


Ec?' 


(14) 


Comparing the strength of the Reynolds stresses and the viscous 
stresses, we can identify Rrp = -(Trp and calculate an effective 
viscous OTv from this equation. This relation is given by 


Q!y 


Rr. 




(15) 


where £ is the turbulent scale height in the BE, which is approx¬ 
imated to be the smaller of the disk scale height H and the radial 
pressure scale height (see |Papaloizou & Stanley] 1986| |Popham| 
|& Narayan|1995| ). 


2.4. Transport of mass and angular momentum 

The mass and AM transport in the disk are two important quan¬ 
tities when we are interested in the effectivity of the Reynolds 
stresses in the BE. Both result from the ID stationary hydrody¬ 
namic equations, where they play the role of eigenvalues. The 
mass flux is obtained from the continuity equation and reads 

M = 2nr'Lur. (16) 


The angular momentum flux follows from the momentum equa¬ 
tion in azimuthal direction and is given by 

j — Inr^Yuriip — Inr^cfrp -I- Inr^Rrp. (17) 

Most commonly, the AM flux J is expressed as a dimensionless 
parameter j, which is the AM flux normalized to the advective 
AM flux at the surface of the star: 


j^j, A = MRiTl^iR.). 


(18) 


Either flux is defined such that if mass or AM are traveling out¬ 
ward, M and j are positive. The three terms in Eq. (17 1 can be 
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identified as the advective AM transport, the viscous transport, 
and the transport due to Reynolds stresses. Both M and j are 
constant with respect to space and time only if the system has 
reached a stationary state, i.e. djdt - 0. Then, M equals the 
imposed mass flux that is set as a boundary condition. Eor non¬ 
stationary states, M and j represent the local mass and AM flux 
at a certain time. 


3. Numerics 


3.1. General remarks 

We used the code FARGO 


with modifications of 


Baruteau ( 2008|l for the simulations presented in this paper. 


FARGO is a specialized hydrodynamics code for simulations of 
an accretion disk in cylindrical coordinates {r,(p) that exploits 
the advantages of the EARGO algorithm ( |Masset|[2000| l. The 
EARGO algorithm is a special method for the azimuthal trans¬ 
port in differentially rotating disks that is able to speed up the 
simulations by up to one order of magnitude and exhibits a 
lower numerical viscosity than the usual transport algorithm. 
The speed-up works best when the velocity perturbations are 
smaller than the bulk flow in azimuthal direction. Since the aver¬ 
age azimuthal velocity Up{r) changes significantly in the course 
of the simulations, the background rotational profile for the 
EARGO algorithm is updated at each time step in order to gain 
full benefit of the method. The algorithm has also been shown to 


produce valid results when shocks are created in the disk (Masset 
|2000)l. The radiative diffusion in the flux-limited approximation 
(Eq. 0 is solved implicitly using a successive over-relaxation 
(SOR) method (see [Muller & Kley 2013| l. The code has been 
tested and used extensively for a variety of applications. In order 
to verify its suitability for simulations of the BE, we compared 
it to the ID BE code described in [Hertfelder et alT]p013] l and 
found a perfect agreement for the ID models. 


3.2. Boundary and Initial conditions 

To model the BE with a net mass flow through the disk we 
have implemented specific boundary conditions in FARGO. At 
the outer boundary we impose a mass flux M which goes into 
the disk and whose absolute value can be chosen as a parameter. 
We require the azimuthal velocity to be Keplerian at the outer 
boundary, M^(rout) = routDK(rout)- Modeling the inner bound¬ 
ary, i.e. the beginning of the star, is more difficult. Since we are 
treating the problem in cylindrical geometry in the disk plane, 
we cannot simultaneously simulate the outer parts of the star. 
Therefore, we start the simulation domain of all our simulations 
at rjn = 1R», where R, is the radius of the chosen star, and set 
the azimuthal velocity to the velocity of the rigidly rotating star, 
Upfr-m) - RtQt, i.e. we impose a no-slip boundary condition in 
i/j-direction. The radial velocity is set to a certain fraction T of 
the Keplerian velocity 

Mr(nn) = -TrinGK(rin). (19) 

This is obviously an outflow BC for the radial velocity at the 
inner boundary and it is necessary since the matter which is in¬ 
serted in the domain from the outer boundary must be allowed 
to leave again at the inner boundary. The factor T determines 
how deep the inner boundary rin is located inside the star: the 
smaller T is, the deeper the domain ranges into the star. The ac¬ 
tual choice of !F is a trade-off between containing a sufficient 
part of the stellar envelope in the domain and avoiding densities 
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and temperatures that are too high at the inner boundary to pre¬ 
vent the simulations from being slowed down too much. A value 
of ^ = 10'® is used in our simulations. We have tested differ¬ 
ent values of T and found that it does not have any influence 
on the initial models apart from shifting the radial prohles by a 
small amount in r. In particular, it does not have an effect on the 
mass flux through the inner boundary once the initial model is in 
an equilibrium state. Instead, the density adjusts such that M at 
the inner boundary equals the value at the outer boundary. Apart 
from the special boundary treatment mentioned above, we im¬ 
pose zero-gradient or Neuman-type boundary conditions for the 
remaining variables, which means that the normal derivative at 
the boundary vanishes. 

The 2D simulations were started from fully relaxed ID ax- 
isymmetric BL models which were generated using the same 
methods as in [Hertfelder et al.| P013] l. Those simulations use 
a classical a-parametrization for the viscosity with a uniform 
value of a = 0.01 both in the disk and the BL. We employed 
these models to have realistic profiles for the density and tem¬ 
peratures as initial conditions. After the ID models reached a 
stationary state we extended them to 2D input data for FARGO and 
ran them again with a small number of azimuthal cells until they 
were fully relaxed. Following this burn-in procedure we started 
the actual simulations from these axisymmetric initial conditions 
with full resolution in r and (f. Prior to the start, the radial veloc¬ 
ity was randomly perturbed by 1% of the Keplerian velocity at 


the surface of the star. The initial profiles are shown in Fig. 10 
below (dashed lines). 


3.3. Model parameters 


In this paper we focus on the non-axisymmetric instabilities and 
the AM transport in the BL. We chose the physical system of a 
BL around a young star, because the ratio of the mass and the 
radius of the star, is much smaller than, for instance, for 

a white dwarf, which we considered in a previous paper about 
the BL. This ratio determines the radial scale height in the BL 
and the smaller H is, the higher the resolution must be in order 
to resolve the BL correctly. Therefore, the computational costs 
increase enormously when one goes to more compact objects. 
However, we expect what we find here to be the generic case. 

The young star we consider in this paper has one solar mass, 
M* = Mq, and is three times the radius of the sun, - 3Rq. 
It has an effective temperature of T* = 3000 K which yields a 
luminosity of 3Lq. The star is not rotating, i.e. = 0, which is 
slightly below the observations ( jCohen & Kuhi||T979| l, but was 
chosen for the sake of a large velocity drop. We impose a mass 
accretion rate of M - 10'®Mo/yr, which is a typical value for T 
Tauri stars such as HL Tau ( Hayashi et al.|199^|Lin et al.|1994| l, 
and is two orders of magnitude smaller than the accretion rate 
FU Orionis stars can reach in outburst phases (Hartmann et al. 
|1993|l. We took a - 0.01 for the viscosity parameter in the disk 


and a = 0 for the BL (see Sect. 2.2 1 . We consider the following 
four simulations for this work. 


- Simulation A: It has a resolution of 512 x 512 and the radial 
domain ranges from to 2/?». This simulation features the 
lowest resolution of all runs performed. 

- Simulation B: This is the reference simulation, and unless 
stated otherwise, we always refer to it in the text. It has a 
simulation domain of [/?,,2f?t] and a resolution of 1024 x 
1024. 

- Simulation C; Here, the number of grid cells is doubled in 
each direction and hence the resolution is given by 2048 x 


2048. The simulation time for this resolution is very high 
and we are not able to study the long-term evolution for this 
run. 

- Simulation D; In order to ensure that the radial domain in 
simulations A-C is not too small, we increased it to [R,,5Rt]. 
While the first three simulations utilized an arithmetic grid, 
we chose a logarithmic spacing in radius for this run along 
with a resolution of 1024 x 1024. 


4. Boundary layer instability 

In this section, we will summarize the key features of the sonic 
instability and its associated acoustic waves in the BL (jGlatzelj 


1988| Belyaev & Rafikov|2012||Belyaev et al.|2012 2013a| l and 
demonstrate that our models are also prone to the sonic instabil¬ 
ity. 


4.1. Sonic instability 

The sonic instability arises in supersonic shear flows and is a 
non-local instability similar to the Papaloizou-Pringle instability 


( Papaloizou & Pringle|1984] Narayan et al.|1987t Glatzel|1988] l. 

It is the analogue of the classical KH instability in the regime of a 
supersonic drop in the velocity of a compressible, stratified fluid. 
The sonic instability has recently been studied in the context of 
the BL extensively by [Belyaev & Rafikov| ( |2012] ). 

There are two ways by which the sonic instability operates: 
The first is the overreflection of sound waves from a critical layer 
that corresponds to the corotation resonance in a rotating fluid, 
i.e. the radius where the pattern frequency of the mode matches 
the rotation frequency of the bulk flow. This mechanism results 
in a dispersion relation that features several sharp peaks and 
hence discrete modes are clearly favored, especially if the den¬ 
sity contrast across the shear interface is high (see e.g. panel (n) 
of Fig. 3 in Belyaev & Raflkov|2012| l. Each mode can be associ¬ 
ated with a pseudo-energy, or conserved action (angular momen¬ 
tum), which can leak through the critical layer when a mode is 
reflected at the corotation resonance ( Narayan et al.|1987 1 . The 
tunneled wave has the opposite sign of the action and, since it is a 
conserved quantity, the reflected wave must undergo an increase 
in its action in order to globally compensate for the additional 
negative action. Thus, the reflected wave is amplified. Further¬ 
more, if a mode is trapped such that it performs multiple reflec¬ 
tions at the corotation resonance, this amplification mechanism, 
called the corotation amplifier (e.g. [Mark[1976[ ), eventually leads 
to instability. 

The other destabilizing mechanism is the radiation of en¬ 
ergy away from the BL. If a mode inside the BL has a negative 
action density, it becomes even more negative through the radi¬ 
ation mechanism and, consequently, is amplified, which results 
in further instability. The dispersion relation of this destabilizing 
mechanism is a smooth function of the wave vector. Therefore, 
a continuum of modes are associated with it. This mechanism 
dominates the overreflection for a small density contrast on ei¬ 
ther side of the interface (see again Fig. 3 of[Belyaev & Raflkov[ 

Each mechanism, the overreflection and the radiation, is 
most efficient in the regime of large density contrast or equal 
density, respectively, and they operate together. The growth rates 
of the sonic instability are approximately l/P*, and the depen¬ 
dence of the growth rate on the density contrast is small, de¬ 
spite the different destabilizing mechanisms. This further dis¬ 
tinguishes the sonic instability from the KH instability. In the 
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Fig. 1. The total radial kinetic energy as a function of the very first 
orbits for three simulations mentioned in Sect. |3.3[ which differ only in 
their spatial resolution. At around orbit 5, the sonic instability starts to 
operate in all three cases, and as a result, the radial kinetic energy rises 
rapidly. The growth rates were obtained by applying an exponential fit 
to the first exponential increase. 


numerical simulations of the BL, a very high numerical res¬ 
olution is necessary to cotTectly reproduce the growth rate of 
the sonic instability. Such a high resolution is computationally 
very demanding and costly, especially for investigations of the 
long-term evolution of the BL. However, as has been shown by 
|Belyaev & Rafiko^ ( |2012| l and revisited in |Belyaev et al.| ( [2012| l, 
a lower resolution merely leads to a slower growth of the insta¬ 
bility. 

We confirm this point by comparing the growth rates of three 
simulations with different spatial resolution. Figure shows the 
evolution of the radial kinetic energy when the sonic instability 
starts to set in for the three resolutions 512 x 512, 1024 x 1024, 
and 2048 x 2048. Indeed, the simulation with the lowest spatial 
resolution underestimates the growth rate by a factor of ~ 1.8. 
However, the simulations with 1024^ and 2048^ cells both yield 
an almost identical growth rate of 1.09 per orbit. The shear layer 
with positive gradient of Q is initially resolved with approxi¬ 
mately 150 and 300 cells in radial direction, respectively. In sim¬ 
ulation A, we only have about 75 cells across the velocity drop. 
We also verify by Fig. [T] that the growth rates of the sonic insta¬ 
bility, which are of the order of 1 per orbit, are actually very high. 
We will now focus again on our reference simulation, although 
all simulations behave basically the same. After about ten orbits, 
the rise of the radial kinetic energy is interrupted briefly and then 
continues with a slightly different growth rate. The paused in¬ 
crease in the instability and the change in growth rate is due to 
the change in density and velocity in the shearing region with 
time. Since the mass inside the BL is delivered from the disk, 
the density in the vicinity of the velocity drop rises. The veloc¬ 
ity itself changes as a result of the initiated angular momentum 
been mentioned and is elucidated in 
, the growth rate of the sonic instabil¬ 
ity shows a weak dependence on the density ratio and the modes 
depend on the slope of the velocity drop, as well. Thus, a change 
in the instability mechanism can be observed at this time. 

We now have a closer look at how the sonic instability man¬ 
ifests itself in the radial velocity and consider Fig.|^for this pur¬ 
pose. The upper panel is a snapshot taken after five orbits and the 
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transport. As has already 


Belyaev & Rafiko^(2012 


lower panel shows Uy ait - 10 orbits. In both panels, the veloc¬ 
ity is normalized to the Kepletian velocity at the surface of the 
star. Figure]^) features a distinct boundary at rcut ~ 1.19, which 
coincides with the inflection point in Q(r) and we can clearly see 
sound waves on either side, propagating away from the bound¬ 
ary. The pattern speed equals the rotation rate of the unperturbed 
flow at the location of the boundary and thus we associate it with 
the critical layer, or corotation resonance, which we have already 
introduced. It is also clearly noticeable that the modes undergo 
a phase shift across the critical layer. Both features are a clear 
signature of the sonic instability and we can compare Fig. [^) to 
Fig. 2 of Belyaev et al.| ( 20121 and Fig. 5a) of Belyaev & Rafikov 


P012| l to verify that our model is indeed unstable to the sonic in¬ 
stability. 

Figure]^) depicts the radial velocity five orbits later and the 
appearance has already noticeably change. The amplitude of the 
sound waves in the lower panel is at least a factor of 2 larger 
than that in the upper panel and amounts to a large percentage 
of the azimuthal velocity. The picture shows a distinct pattern of 
wave front crossings and reflections for r < rcrit and the critical 
layer at rmt = L19 has already begun to smear out at this point 
in time. When the color scale is saturated we recognize three 
regions where the amplitude of the waves is dropping owing to 
the wave fronts crossing each other. They lie at r « 1.16,1.12, 
and 1.07. Inside the corotation resonance, the waves are reflected 
near the inner boundary and at the critical layer, and are thus 
trapped between the two radii. We have traced three wave fronts 
with the black lines to illustrate the pattern that is exactly the 
overreflection mechanism that we have already elucidated. The 
pattern number of the sonic instability (15 for the simulation 
shown in Fig. does not depend on details of the inner bound¬ 
ary condition. This has been verified by running the same simu¬ 
lation with identical perturbations with modified boundary con¬ 
ditions, where all variables are kept at their initial values for all 
time (“do-nothing” BCs). The pattern number is, instead, set by 
the initial random perturbation. For a different seed of the ini¬ 
tial perturbations the pattern number will be different. When the 
trapped waves undergo a reflection at the critical layer, part of the 
wave action leaks through the critical layer and has the opposite 
sign compared to the reflected wave, whose amplitude grows as 
a consequence of the global conservation of wave action. Thus, 
the whole process is self-sustaining and even grows with time. 
The leakage of wave action is also clearly visible in Figure [^) 
for r > Tcrit. Now the reason for the change in the growth rate of 
the sonic instability after about ten orbits (see Fig. becomes 
apparent; It is simply the overreflection mechanism setting in, 
gaining in strength, and growing dominant. 


4.2. The BL modes 

In Section [4~T| we described the linear growth regime of the sonic 
instability and analyzed its manifestation in our models. After 
only 20 - 30 orbits the sonic instability saturates and the sys¬ 
tem remains in the excited state with a radial kinetic energy that 
is approximately three orders of magnitude larger than before 
the sonic instability set in. There are still considerable variations 
in the radial kinetic energy (see Fig. [^; however, we suggest 
that they a result of a secondary KH like instability that sets in 
periodically and dies out shortly after. This theory is examined 
thoroughly in Section [5^ 

After the sonic instability has saturated we observe waves 
that are acoustic in nature propagate through the largest part of 
the domain. Those modes have been excited by the sonic insta¬ 
bility and do not die out for the whole simulation time, which 
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Fig. 2. r-(f plot of the radial velocity in units of uy,{R,) at f = 5 (a) and f = 10 (b) orbits. The velocity in the upper panel, a), has been multiplied 
by a factor of 2 in order to use the same colormap as below. Both snapshots show the linear stage of the sonic instability and the middle branch in 
action. The sound waves are propagating toward the inner and the outer edge of the domain in an oblique fashion. The three black lines in panel 
b) illustrate the wavefronts. The arrows denote the direction of the unperturbed azimuthal flow. 


amounts to over 2000 orbits. Considerable effort has been ex¬ 
pended by [Belyaev & Rafikov] ( |2012| l and [Belyaev et al.| ( [2012[ 
2013a[ l in order to analyze these acoustic modes by means of 
linear stability analysis and derive proper dispersion relations. 
Since this undertaking is quite involved we limit ourselves to re¬ 
viewing those elements that are crucial for the understanding of 
our results. 

One finds that the acoustic modes of the BL are related to the 
waves in a plane parallel vortex sheet in the supersonic regime 
( Belyaev et al.[[^13a[ l. For a discontinuous drop of velocity in 
the isothermal vortex sheet in absence of stratification or rota¬ 
tional effects it can be shown that the linearized Euler equations 
yield three different expressions for the dispersion relation u, 
which are denoted the lower, middle, and upper branch ( [Belyaev[ 
[& Rafiko^[2012| l. The three branches differ in the direction of 
propagation of the modes on either side of the interface (ve¬ 
locity drop) and their general appearance is depicted in Fig. 4 
of [Belyaev et al.[ ( [2013a[ l. The lower and upper branch are very 
similar to each other. They both show a region where the waves 
propagate parallel to the interface and a region where the wave 
vector forms an angle between 0 and n/l with respect to the nor¬ 
mal vector of the interface. Transferred to the formalism of our 
models, the modes of the lower branch propagate purely in (p- 
direction in the BL and both in r and in ip in the disk, where 
krlk,p » 1, however. For the upper branch, the opposite ap¬ 
plies. Modes of the middle branch propagate in an oblique man¬ 


ner both in the disk and in the BL with the same wave vector. 
Furthermore, the waves undergo a phase shift of n at the inter¬ 
face. This holds true for all branches. The dispersion relations of 
the three branches of the vortex sheet can then be extended by 
including radial stratification and the effects of rotation (Corio¬ 
lis force) and one can find approximate dispersion relations for 
a more realistic case of a BL. This was done in [Belyaev et al.[ 
P013a[ ) for the isothermal BL. It would be interesting to see how 
a realistic equation of state and the radiation transport affect the 
dispersion relations. However, this is beyond the scope of this 
paper. 

Figure]^ is a r-tp plot of the radial velocity after 45 orbits. 
Three white dashed contour lines of the azimuthal velocity Uy, 
have been overdrawn. All velocities are, as usual, given in units 
of the Keplerian velocity at the surface of the star. The image 
shows a standing wave in azimuthal direction for radii < 1.17 
and a wave that propagates into the disk with wave vector com¬ 
ponents in both r- and i/j-directions. Additionally, the amplitude 
of the wave changes sign at r 1.17, i.e. the phase is shifted by 
7r/2. This pattern is a clear signature of the lower branch of the 
acoustic modes, which we can observe in action here. The con¬ 
tour lines of reveal that the interface between the disk and the 
BL, or in other words the region where the velocity starts to de¬ 
viate substantially from Keplerian, is strongly deformed. The de¬ 
formation of the interface is evidently due to the large scale vor¬ 
tices that reside at the base of the BL (see Fig.[^. At r = 1.2 the 
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Fig. 3. Two-dimensional plot of the radial velocity at t = 45 orbits. 
The direction of the unperturbed flow is from bottom to top, i.e. in the 
direction of growing ip. The dashed white lines are contour lines of the 
azimuthal velocity at the levels 0.15,0.3, and 0.7 (from left to right). 
All velocities are given in units of uy^iR,). The radial velocity develops 
a pattern with vortices in the BL and shocks at the top of the BL. This 
pattern most likely resembles the lower branch of the acoustic modes 
discussed in|Belyaev et al.|([2013a|>. 


sound speed amounts to ~ 8% of the Keplerian velocity, whereas 
the gas rotates with 70%MK(f^*)- Therefore, the azimuthal mo¬ 
tion is still nearly ten times super sonic and information about 
the interface deformation cannot propagate upstream. Hence the 
oblique shocks, which are clearly visible in the radial velocity, 
have developed. The whole system described above rotates in 
a prograde way with a pattern speed of Qp = ot/m, where w 
is the oscillation frequency in the inertial frame and m the az¬ 
imuthal wavenumber. For the episode described here (orbit 45), 
Qp ^ 0.3. 

The waves of the lower branch, which we find almost exclu¬ 
sively in our simulations, are tightly wound with a leading spiral 
and we expect kr/m » 1, where k and m are the radial and az¬ 
imuthal wave numbers, respectively. Therefore, one can apply 
the WKB approximation and assume perturbations of the form 

6f oc exp [/ {kr + rmp - wf)]. (20) 


A linear stability analysis for a fluid disk with a polytropic equa¬ 
tion of state then yields the dispersion relation 


2„2 


{u — mQ) = ic + k c\ 


( 21 ) 


where k - ^(r^Q)) is the epicyclic frequency and atlm - 

Qp (|Binney & Tremaine||2008|. Of course, Eq. (pT]i is only an 
approximation to the dispersion relation of the waves in our 
fully radiative disk. However, it will mainly be used to explain 
the wave dynamics in the disk. Furthermore, the WKB approx¬ 
imation is no longer valid in the BL where we have seen that 
kr ai 0. In the disk we can assume a Keplerian rotation pro¬ 


file, Q 


-1.5 


and consequently the epicyclic frequency is 


real. Therefore, there is a region around the corotation resonance 
Q = Qp, flanked by two Lindblad resonances Q = Qp + Kim, 
where k is imaginary and the waves are evanescent. Formally, 
according to Eq. (|2T]l the wave number k becomes zero at the 
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Lindblad resonances. However, the WKB approximation breaks 
down at these points since the assumption krjm » 1 is no longer 
valid. Since the whole pattern of waves and interface deforma¬ 
tions rotates with a pattern speed 0 < Qy^ < 1, there are two 
corotation resonances, one in the BL and one in the disk, ei¬ 
ther of them flanked by two Lindblad resonances. Assuming 
the disk rotates with Keplerian frequency, we can easily derive 
the corotation and Lindblad resonance, which are then located 
at Tcor = Qp^'^^ and rund = (Qpm/(OT + 1))“^'^^. The waves or 
the weak shocks we discussed earlier launched at the BL now 
propagate into the disk until they reach the Lindblad resonances, 
where they are reflected and propagate back toward the BL. In 
the BL the waves might be reflected by Lindblad and corotation 
resonances as well; however, Eq. (j2T| is not applicable there. 
Through this process of consecutive reflections the waves can 
become trapped between the BL and the disk. 

This idea of the trapped modes has been proposed by 
[Belyaev et al.| (|2012| l and was indeed observed in their simu¬ 
lations. In Fig.j^we do not And such clear evidence for reflected 
shocks. Reflection of the wave depends, however, on whether 
the modes are absorbed before they reach the evanescent region 
in the disk. We observe that the modes in our simulation can 
only travel until r 1.5 before they are substantially damped 
by shocks, radiation, or viscous damping due to the a-viscosity 
applied in the disk (see Sect. |2.2| i. We have run simulations with 
larger domains in order to investigate a possible influence of the 
outer boundary on the reflection of the waves. These test runs 
confirm the early damping of the waves observed in the simula¬ 
tions presented in this work. The effective viscosity in the disk is 
likely to play a major role in this damping process. The actual in¬ 
teraction between the waves and the MRI turbulence in the disk 
can, however, only be investigated in real MHD simulations. 

In Fig.|^we plot the velocities Ur and as a function of (r, (p) 
at f = 173 orbits, a more energetic state than the one depicted in 
Fig.[^(see Fig.|^below). It is immediately apparent that the state 
at f = 173 orbits is far more violent than the situation shown in 
Fig. [3 The amplitude of the radial velocity has increased by a 
factor of 2 and the shocks have grown more intense and vigor¬ 
ous. While the general pattern of the radial velocity remains un¬ 
changed and still constitutes the lower branch with the vortices 
at the base of the BL and the outgoing waves at the top of the BL, 
the wavefronts are now heavily deformed by the shocks. These 
shocks are a consequence of the strongly distorted interface as 
can be seen in Fig.|^). The lines of constant u^p in the r-tp space 
(the contours) have a wave-like shape with many irregularities 
and the transition from zero rotation rate to Keplerian rotation 
is spread over ~ 0 . 2 /?* and has broadened considerably. Though 
the radial velocity is more chaotic in Fig.|^) than in Fig.[^ there 
is still a dominant mode visible and the azimuthal wave numbers 
are m = 14 and 24 , respectively. As we shall see later, the an¬ 
gular momentum transport by the waves is very efficient in this 
stage. 

For completeness, in Fig. we show a r-i,p-plot of the radial 
velocity that spans the whole simulation domain. The snapshot 
was taken after 200 orbits while the system was in a quiet state. 
The pattern of the lower branch with its vortices at the base of the 
BL is clearly visible. We can also see that the waves emerging 
from the top of the BL are indeed tightly wound and that the 
spirals are leading. Figure|^shows another azimuthal wave with 
m -3 that we had not detected before, which is still a little faint 
is this figure. Interestingly, [Belyaev et al.| ( |2012| l also And such a 
large-scale m - 3 global mode in their full In simulations. The 
wavenumber of the high-frequency pattern at the base of the BL 
is given by m - 24 
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Fig. 4. r-if plot of the radial velocity (a) and the azimuthal velocity (b) at t = 173 orbits. Both quantities are given in units of Uy^R,). The dashed 
white lines in both plots trace the contour of = 0.15 and 0.6 (from left to right). The direction of the unperturbed flow is from bottom to top. The 
images have been taken at a high activity state and the strong shocks in and the heavy deformation of the interface (see b) are clearly visible. 
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Fig. 5. Full-disk plot of att = 200 orbits. The direction of the unper¬ 
turbed flow is counterclockwise. The radial velocity is given in units of 
Uy(R,). 
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Fig. 6. Temporal evolution of the radial kinetic energy of our reference 
simulation. The time is given in units of the Keplerian orbit at r = R„ 
P,. The switching between quiet and outburst states described in the 
text is clearly noticeable. We have marked the exact times for some of 
the outbursts for easier comparison with the text. 
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The middle branch is a transient mode and dominates during 
the linear growth of the sonic instability. It can be observed in 
the early stages of the simulation and the oblique wavefronts on 
either side of the interface are clearly visible in Fig.|^ Both panel 
a) and b) show features of the middle branch, yet they are more 
distinct in b). The upper branch, on the other hand, might be 
observed during the late times of the simulations. We defer the 
discussion of this issue to Sect. 1531 


5. The effects of the instabilities 

In the previous section we discuss the hydrodynamical instabil¬ 
ities in the BL in detail. Now we show the implications of the 


BL instabilities and the physical results of our simulations of the 
radiative BL around a young star. 

5.1. Temporal evolution of the BL 

Figure visualizes the total radial kinetic energy, Fkin.iad = 
for which we summed the radial kinetic energy of 
every cell in the simulation domain. For a given time, the radial 
kinetic energy reflects the strength of the instabilities. There¬ 
fore, we can monitor the time-dependence of the instabilities 
by tracking Fkin.rad- The reference simulation shows a strongly 
time-dependent behavior that can be characterized as a recurring 
series of states of outburst and quiescence. When the kinetic en- 
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ergy reaches a maximum, the instabilities are vigorous and the 
interface between the BL and the disk is heavily deformed (see 
Fig.g. In the quiescence states, on the other hand, the amplitude 
of the velocity perturbations drops considerably, often by a fac¬ 
tor of more than 10, and the interface deformation is less strongly 
pronounced. The strength of the outbursts does not follow a clear 
pattern. In the beginning, the outbursts get stronger, with the sec¬ 
ond outburst at f ^ 173 orbits being the most powerful. During 
this outburst the amplitude of the perturbations in the radial ve¬ 
locity reaches over 20% of the Keplerian velocity at the surface 
of the star. Thus, the amplitude of the radial velocity becomes 
comparable to that in azimuthal direction, or even exceeds it for 
r < 1.1. The third outburst at f ^ 280, which is less intense than 
the second one, is followed by a longer period of less activity 
or quiescence during which some small peaks suggest the occur¬ 
rence of very weak micro-outbursts. After more than 500 orbits, 
the system starts to go through a series of outbursts again. 

The high and low instability states mentioned above corre¬ 
spond to episodic phases of high and low activity in mass accre¬ 
tion and angular momentum transport, as can be seen by com¬ 
paring Figs. 1^ and Figure 7 depicts the time-dependence of 
the mass accretion rate, Eq. and of the angular momentum 
transport, Eq. ( [T7| ), for all radii inside the computational domain. 
In phases of outbursts, which are marked by the vertical dashed 
lines, the mass flux through the disk is higher than in the qui¬ 
escence state. In certain regions, M increases by more than two 
orders of magnitude. In the high activity phases, M distinctly 
rises in the outer parts of the simulation domain as well. The an¬ 
gular momentum transport, j, shows a similar trend. At times of 
strong instability, the transport of angular momentum, which is 
mainly directed to the center of the disk, also grows remarkably, 
quite often by a factor of about 100. 

Unlike in a stationary state where mass and AM are trans¬ 
ported to the center of the disk (see Eq. [T7] i, we observe that 
the transport of mass does not have a preferred direction in cer¬ 
tain regions of the simulation domain. We consider for this pur¬ 
pose the upper panel of Eig. [7] There is a clear boundary at 
r 1.2 - l3Rt, which is approximately given by the location 
of the maximum of Q and divides the rather organized, inward 
directed mass flux in the disk from the alternating transport in 
the BL, where mass is sloshing about back and forth in radius 
because of the passage of waves which have alternating signs 
of mass transport. This is clearly visible in the form of the di¬ 
agonal features in the upper panel of Eig. |7] Thus, the direc¬ 
tion of the mass transport in the BL is both a function of time 
and radius, and the variability with time is on the order of one 
orbit. After roughly 600 orbits, the strength of the mass trans¬ 
port ceases slightly, only gaining strength again during the sub¬ 
sequent outburst periods. This development matches that of the 
radial kinetic energy (Eig.|^ where we can also detect a change 
in activity starting after 600 orbits. After 1200 orbits another ma¬ 
jor change in the mass transport pattern takes place. The whole 
system seems to calm down and the sloshing behavior with its 
frequent changes in direction ceases. 

The AM flux does not feature a variability as intense as the 
mass flux; however, the diagonal features that have also been 
discussed in the context of M are still present. An analysis of 
the individual components of j confirms that these features arise 
from the advective part of the AM flux, which shows a sloshing 
behavior similar to the mass flux. Apart from these features, j 
assumes mostly negative values (see the lower panel in Eig. |^, 
meaning that AM is transported inward. This evidence agrees 
with our expectations that in the disk j is negative for the rea¬ 
sons pointed out earlier. One very distinct feature in the graph 
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Fig. 8. The velocity in azimuthal direction in units of the Keplerian ve¬ 
locity at the stellar surface for seven different orbits. The curves demon¬ 
strate the spin-up of the gas for r < 1.2 and the spin-down for r > 1.2 
over time. 


of the AM flux is the very pronounced boundary that starts at 
the beginning of the simulation at r close to \.2R, and moves 
inward up to r 1.025/?* during the first 800 orbits. This line 
can be characterized by an enhanced AM flux in its near vicin¬ 
ity and a shift in direction from inward (r < rune) to outward 
(f > nine)- Eigure]^ visualizes the destination of the AM which 
is accumulated at this interface as a result of transport from both 
sides. It is spent to spin up the gas at r < 1.2. Hence, the BL is 
broadened considerably during the course of the simulation. The 
AM needed for the spin-up of the inner layers is extracted from 
the faster rotating part at r > 1.2, which is, consequently, slowed 
down. By comparing Fig. and Fig. we can deduce that the 
AM transport is the strongest where the gradient of Q. is highest, 
i.e. where the shearing is strongest. 

It can be recognized from Fig. that AM is mostly trans¬ 
ported inward although the mass flux is also directed outward. 
Especially when M > 0 and large, one would expect a consid¬ 
erable amount of AM to be transported outward along with the 
matter. Thus, one possibility for the AM flux to still remain neg¬ 
ative is that the Reynolds part of Eq. 0 is negative and larger 
than the advective part (the viscous part is zero in the BL). How¬ 
ever, an analysis of several time steps showed that very often 
the advective AM transport is negative and large while the mass 
flux is positive and large. This apparent discrepancy between the 
mass and the AM flux is due to an anti-correlation between the 
radial and azimuthal velocity. The mass accretion rate periodi¬ 
cally points in- and outward, depending on azimuth, owing to 
the vortices of the radial velocity. The sum of M over the whole 
azimuthal domain for a given radius can therefore be positive 
(meaning that mass is transported outward). However, accord¬ 
ing to Eq. 0’ the total advective AM flux can easily remain 
negative if is smaller at azimuthal locations of positive M 
than where it is negative. This is indeed very often the case and 
we have picked an illustrative example to stress this point (see 
Fig.|g. It is a snapshot of the mass accretion rate, the radial ve¬ 
locity, and the angular velocity at orbit 145, where we plotted 
the dependence on azimuth for r = 1.12/?,. The figure clearly 
shows that whenever M{ip) is positive and large, u^p{(f) is small. 
The comparison between the radial and the azimuthal velocity 
reveals that the anti-correlation between Ur and is such that 
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Fig. 7. The upper panel shows the temporal evolution of the mass accretion rate in the radial direction, M{r, t), in units of solar masses per year. It 
is obtained from the mass accretion rate per cell (see Eq.[T^ by integrating over the azimuthal direction. The lower panel likewise visualizes the 
AM flux normalized to the advective AM flux at the stellar surface (see Eq.|18|l, which is also integrated over azimuth. We note that we applied 
a so-called symlog scaling for the plot, which means that we have a logarithmic scaling for both the positive and the negative values. Only in the 
vicinity of 0 is the scale linear. A negative value means that mass or AM is transported to the center of the disk. The phases of strong activity in 
the radial velocity (see Eig.[^ marked with the black dashed lines here) correspond to phases of increased mass and AM transport. The diagonal 
features in the upper panel are a clear signature of waves, which have alternating signs of mass transport. This sloshing behavior spreads to the 
AM flux (lower panel) via the advective component of j. 


an extreme of the azimuthal velocity coincides with a root of the 
radial one. Therefore, the two patterns that are similar, although 
not perfectly equal, are shifted by one quarter of the azimuthal 
wavelength of the m — 1 mode at this annulus. 


5.2. Outburst sequence 

We now take a closer look at the cause of the periodic outbursts 
and the behavior of the primitive variables in the disk during 
such an outburst and we consider Fig. [^for this purpose. The 
graphic shows the velocity in azimuthal direction, the surface 
density, the pressure, and the temperature in the midplane of the 
disk for five different times that are centered around the outburst 
at f = 173 orbits and are not equidistant. These time steps cover 
the second outburst in Fig. This part of the radial kinetic en¬ 
ergy evolution is shown as a zoom-in in the lower right panel of 

Fig.[Tg 

At the beginning of the outburst, the surface density has in¬ 
creased quite strongly, especially in the vicinity of r « 1.25 
where it is almost two orders of magnitude higher than at the 
beginning of the simulation. This pile-up of mass is a result of 
the fact that mass is constantly entering the BL from the outer 


disk, but the instabilities in the BL are not yet efficient enough 
to guide it all through. This situation is universal for all ob¬ 
served outbursts, which are preceded by a phase of quiescence 
during which the activity and the mass transport in the BL are 
not as strong as in outburst phases. Therefore, the whole pro¬ 
cess described here using the example of the first outburst can be 
equally applied to the other outbursts. During the accumulation 
of the mass just outside the BL, the azimuthal velocity profile is 
smooth and quasi-stationary. The high amount of mass leads to 
a high temperature and hence a high pressure, which means that 
Ug, is remarkably sub-Keplerian in the first half of the domain 
since the gas is stabilized by pressure. 


After a sufficient amount of mass has accumulated, a run¬ 
away process sets in during which and the other quantities 
change very rapidly. According to Figs. 7]and[T^a great amount 
of mass is guided through the BL onto the star over such an out¬ 
burst. The process is initiated outside the BL at about r « 1.3. 
The gas is both accelerated and decelerated in azimuthal direc¬ 
tion by the Reynolds stresses in and near the BL, with the bound¬ 
ary between the spin-up and spin-down region moving inward 
with time. Since the effect of the deceleration is much stronger, 
mass moves inward. As a consequence of the decreasing density. 
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Fig. 10. Plots of the azimuthal velocity, the surface density, the pressure, and the midplane temperature for five time steps centered around the 
second outburst. The dashed lines denote the initial conditions of the simulations. The dash-dotted line in the upper left panel gives the azimuthal 
velocity for pure Keplerian rotation. In the lower right panel we have included a zoom-in of Fig.|^for better identification of the exact location of 
the plot times. 


the pressure decreases and the gas is even less stabilized and can 
accelerate its inward movement. After the density is satisfacto¬ 
rily depleted, the Reynolds stresses decrease and the violent ac¬ 
cretion process comes to an end. Then the density grows again 
because of the lower stresses and the azimuthal velocity is slowly 
regains its pre-outburst profile. 

The variation of the midplane temperature during the out¬ 
burst phase is consistent with what we have already discussed 
for the other physical quantities. At the beginning, despite the 
decreasing density, the temperature rises for r < 1.6R» and even 
develops a small bump around r = 1.25/?*. The increase in tem¬ 
perature is due to the dissipation of energy through the growing 
Reynolds stresses and indicates the region where the energy re¬ 
lease by the waves is largest. This radius is in agreement with 
what we expected from the plot of the azimuthal velocity if we 
also take into account that the energy is redistributed by radia¬ 
tive diffusion. Subsequently, the surface density drops more and 
more and the Reynolds stresses decrease as well. Accordingly, 
the temperature also drops, and even drops below the level of 
the initial state. After enough mass has again piled up around the 
BL, the temperature will rise to its pre-outburst level. 

The outburst sequence shown here for the second outburst of 
the reference simulation is exemplary for the general behavior 
of outbursts. What we have just discussed might be less clearly 
visible or even more distinct, depending on the strength of the 
outburst. However, the mechanism is the same throughout. From 
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simulation C, which has a resolution of 2048 x 2048, we find 
that the outburst is far more violent than the one we discussed 
here. The azimuthal velocity frequently develops an even more 
pronounced local maximum and minimum for r < 1.2R, and 
the depletion of mass is more effective. Thus, it is possible that 
the instabilities become stronger with increasing resolution. The 
increasing strength of the outbursts might instead be due to the 
other modes becoming excited through the initial perturbations, 
as is typical for classical KH instability. 

We propose the following scenario for the physical cause of 
an outburst such as the one described above. The radial veloc¬ 
ity has developed weak shocks at the top of the BL. Every time 
a fluid element passes through that shock it is slowed down by 
these oblique shocks. The cumulative effect of all these shock 
passages leads to a decrease in the azimuthal velocity of the gas 
and creates a wide plateau in the velocity profile (see e.g. orbit 
160 in Fig. [T0) i. Such a flat region with a weak velocity gradi¬ 
ent is prone to the classical Kelvin-Helmholtz (KH) instability, 
which sets in and induces the rapid changes in the BL. In ab¬ 
sence of rotation, it follows from Rayleigh’s stability equation 
that a necessary condition for instability is the occurrence of an 
inflection point in the velocity. A stronger form of this condi¬ 
tion was given by |Fj0rtoft| ( | 1 950| l, who showed that the velocity 
profile not only needs to have an inflection point, but also has to 
satisfy the condition u"{u-u^) < 0 somewhere in the flow, where 
Ms is the velocity at the inflection point. However, this require- 
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Fig. 9. The upper panel depicts the mass accretion rate, the middle panel 
the radial velocity, and the lower panel the azimuthal velocity. All three 
quantities are plotted as functions of azimuth at the fixed radius r = 
\.\2R, after 145 orbits. The upper two curves are shifted in phase with 
respect to the lower one such that maxima in M (nearly) coincide with 
minima in u^. 


ment is only true for parallel flows. A normal-mode analysis of 
the linearized Euler-equations of a rotating incompressible flow 
with zero mean radial velocity eventually yields 


{lo + imD.) 


^ 1_5 

dl-2 f Qy 



— im 


1 dZ 
~r^ 


0 = 0 


( 22 ) 


for 2D disturbances. Here, 0(r) is the amplitude of a stream 
function 0'(r, ip, t) - (p{r) exp(wf -i- imtp) such that m' = iimplr 
and = -dip!dr, where the prime denotes perturbed quantities; 
Z = \!r d{ru,^)!dr is the vorticity. Equation (22i is the equiva¬ 
lent of Rayleigh’s equation for a rotating fluidiT&yleigh showed 
that if 0 = 0 at two boundaries R\ and R 2 enclosing the flow, a 
necessary condition for instability is that the gradient of Z must 
change sign at least once in the interval R\ < r < R 2 (|Drazin & 


|Reid|2004| l. We take this analogue of Rayleigh’s inflection point 
theorem as a motivation to study the vorticity and its derivative, 
though Eq. (22 1 is clearly an idealization of the situation at hand. 
However, a detailed stability analysis of the RHD equations of a 
rotating stratified fluid is beyond the scope of this work. 

Eigure [redisplays the derivative of the vorticity, dZjdr, for 
three points in time. At orbit 50, the system is in a state of quies¬ 
cence (see also Eig.l^l. After 150 orbits the second outburst has 
just started to evolve, and 50 orbits later the system is recover¬ 
ing from the outburst. While the vorticity derivative of the quiet 
and post-outburst state does not change sign for r > 1.2R», it 
does change several times for orbit 150. The region r > 1.2/?* is 
of interest here because this is where the rotation profile is flat. 
There are, of course, multiple changes of sign for r < 1.2R, in 
all curves. However, those points are associated with the strong 



Fig. 11. Derivative of the vorticity, dZjdr, as a function of radius for 
orbit 50 (quiet state), orbit 150 (pre-outburst state), and orbit 200 (post¬ 
outburst state). The inset is a zoom-in of the region 1.2 < r < 1.4 where 
multiple changes of sign are clearly visible for orbit 150. 


shearing regions that we have already discussed in Sect. [4 After 
analyzing several time steps before, during, and after outbursts, 
we draw the following picture; The consecutive shock passages 
near the maximum of Q. cause the flow field to adjust such that 
multiple changes of sign are created in the derivative of the vor¬ 
ticity. Since this is a necessary condition for the instability of a 
rotating, hydrodynamical flow, it is probable that a classical KH 
instability develops in the flat regions of the flow. The role, how¬ 
ever, that the KH instability plays in the context of the enhanced 
accretion and AM transport is not yet entirely clear. It might, on 
the one hand, act as a trigger for a mode of the sonic instability 
which is excited during the outburst and efficiently drives the ac¬ 
cretion. Once a sufficient amount of the mass in the BE has been 
accreted, the mode can no longer be sustained and the system 
returns to quiescence. One candidate for this mode might be the 
upper branch, which we normally do not observe in our simu¬ 
lations. A scenario similar to this option was also discussed in 
[Belyaev et al.| ( |2013a[ ). If, on the other hand, the KH instability 
itself is responsible for the enhanced accretion, AM should be 
transported mainly by turbulent stresses, as opposed to the wave 
transport through the mode mentioned above. We did not find 
clear signs of turbulence during the outburst, nor did we find ev¬ 
idence of the upper mode. It is, however, plausible that the AM 
transport is still accomplished by waves and the KH instability 
might trigger or enhance a mode of the lower branch. 

We did find evidence for the existence of a KH instability 
during the outburst. We consider for this purpose Eig. [T^ a 2D 
image of the z-component of the vorticity Vxm at orbit 172.6, i.e. 
at the peak of the second outburst. It shows patterns that are com¬ 
monly associated with the classical KH instability, for instance, 
the curly arms and the cat’s eyes. By analyzing simulation C we 
find that these patterns become more distinct with increasing res¬ 
olution (see also Eig. 11 in Belyaev et al.|2012) l. The KH features 
are, however, masked to some extent by the waves generated by 
the main instability. This becomes clear when we investigate the 
outbursts at later times when the main instability is decreasing. 
In addition, the flat region where we expect KH instability to 
work occurs in the fast rotating part of the disk. Static grid codes 
suffer from disadvantages when the fluid is boosted and fluid 
mixing instabilities can be suppressed (|Springel|2010|l. 
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Fig. 12. Image of the z-component of the vorticity, V x «, at the peak of 
the second outburst at t = 172.6 orbits. The patterns in this image might 
be indicative of a KH-like instability. 



Radius r [i?,] 

Fig. 13. The azimuthally averaged effective temperature of the BL for 
six time steps centered around the second outburst (see Fig. 1^. Teff in¬ 
creases considerably during the outburst and then declines to approxi¬ 
mately the pre-outburst level ~ 30 orbits after the peak of the outburst. 
The black dashed line corresponds to the surface temperature of the 
standard solution of|Shakura & Sunyaev|(|197'3t. 


5.3. Effective temperature 


We now investigate the outburst described in Sect. |5.2| from an 
observer’s point of view and look at the radiation emitted by the 
BL and the disk. In Fig. 10 we have exclusively considered phys¬ 
ical quantities, which are measured in the equatorial plane of the 
disk, among them the midplane temperature T^. If the medium is 
optically thick, however, the emergent spectrum is dominated by 
the temperature on the surface of the disk where the vertical op¬ 
tical depth is approximately one, and locally resembles the spec¬ 
trum of a blackbody radiating with temperature Teff. Since we 
have no knowledge of the vertical structure of the disk, we use 
the approximation of Hubeny ( 1990[), which is a gener alization 
of the gray atmosphere (e.g. |Rybicki & Lightman|1986| l. The ef¬ 
fective temperature can then be calculated with an approximated 
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optical depth in vertical direction ( |Suleymanov|1992| ): 

^eff = T'.VrefF 

3 V3 1 

Teff - o’^R 7-• 


(23) 


Here, tr and tr (t = KpH = O.SkY) are the Rosseland and the 
Planck mean optical depth, respectively. 

Figure 13 shows the effective temperature of simulation B 
calculated according to Eq. ( [23] l. At the beginning of the outburst 
(orbit 160) the effective temperature features only a slight peak 
in the BL and the temperature in the BL and in the disk differ at 
most by 500 Kelvin. The reason why the effective temperature of 
the BL and the disk are very similar at the beginning of the out¬ 
burst episode is the following. The system, which is just on the 
verge of an outburst at this point in time, comes out of a preced¬ 
ing state of quiescence. Therefore, the energy dissipation of the 
waves is small and the density is rather high owing to the ineffi¬ 
cient mass transport. Consequently, both the midplane tempera¬ 
ture and the optical depth in the BL are comparable to the disk 
and Teff is only slightly larger in the BL. Over the course of the 
next 13 orbits, however, the system goes through an episode of 
effective mass and AM transport, and the midplane temperature 
increases and the surface density decreases (see Fig. [T0| . Thus, 
according to Eq. (|23]), Teff increases rapidly in the BL yielding 
a maximum of 6000 K at r « 1.3, which is twice the temper¬ 
ature in the disk. Care should be taken when comparing in 
our models with conventional BL models where an a-viscosity 
has been applied. The latter typically feature a BL that is con¬ 
siderably smaller in radial extent than in the models presented 
here. The luminosity of the BL, however, does only depend on 
the physical parameters of the system and must be the same in 
both models (jGMM//?* for a non-rotating star). Since the lumi¬ 
nosity is the integral of the flux F — over the area of the 

BL, it follows that the thinner the BL, the higher the effective 
temperature must be, and vice versa. The relation between the 
maximum effective temperature in the BL 0 bl and in the disk, 
can roughly be estimated by (|Lynden-Bell & Pringle] 19*^ 


^BL / 6 

= 1.56 -- 
9d \R* 


-1/4 


(24) 


where 6 is the width of the thermal BL (see e.g.|Regev & Bertout| 
|1995[ [Popham & Narayan||1995| l. Equation p4[) yields a factor 
of ~ 2 for (5 = 0.3Rt, which is in good agreement with what 
we observe in our simulation (see Eig. [H- Also visualized in 
Ei g. [T3] (black dashed l ine) is the surface temperature according 
to Shakura & Sunyaev (19731, which is given by 


Teff(r) 


_ I3GM,M 
\87rrVsB 






1/4 


(25) 


where j is the normalized AM flux. We took 7 = 1.1 for Eig. 
which is the mean normalized AM flux at orbit 160 measured 
from the simulation. The effective temperature derived from our 
simulation agrees very well with Eq. ( |25] l in the disk. In the BL, 
of course, the analytically derived formula fails to capture the 
rise of the effective temperature. 

Another important point for the observational appearance of 
the BL is the azimuthal dependence of the effective tempera¬ 
ture, which is visualized in Eig.[^for orbit number 173. We can 
see from this depiction that Teff reaches temperatures of nearly 
10000 K in the simulation domain which is high, considering 
the case of a young star. There is strong variation of T^g in ip- 
direction, especially for the annulus r ^ 1.3, which is due to 
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Fig. 14. 2D plot of the effective temperature of the BL at orbit 173, 
which marks the peak of the second outburst (see Fig. Strongly 

varies in azimuthal direction and reaches temperatures up to 10000 K. 
The curly arms visible for r < 1.2 are indicative of a KH-like instability. 


Fig. 16. The width of the BL as a function of time for the four simula¬ 
tions mentioned in Sect. 13.31 The width is calculated from the radii of 
the region where 0.02 < u^juyi < 0.8. All four simulation settle down 
to a BL width of approximately 20 percent of the stellar radius. Simu¬ 
lation A approaches this value more slowly and in a slightly different 
manner, which is due to the very low resolution of this particular simu¬ 
lation (512x512). 


the shocks of the radial velocity that create density perturbations 
in the surface density and then affect the vertical optical depth 
and the effective temperature. The i^p-dependent pattern of Teff 
rotates with the pattern speed of the acoustic modes, which is 
close to the Keplerian velocity, and will express itself in time- 
dependent variations of the light curve. This is an important 
point because many systems show light curve oscillations with 
periods that are close to the orbital period at the surface of the 
star. In cataclysmic variables, for instance, one observes rapid 
oscillations called dwarf novae oscillations (DNOs, [Warner & 


Robinson]|1972j l and quasi-periodic oscillations (QPOs, Patter¬ 
son et al.||l977^. Both phenomena are not yet fully understood 


and it is tempting to see a connection to the BL modes on the 
basis of the protostar simulations presented here. However, fur¬ 
ther research is necessary to clarify this point. We also note that 
Fig. [I4] shows clear indications of KH instabilities. These fea¬ 
tures are even more pronounced in snapshots of taken before 
the peak of the second outburst and reinforce the point that a 
KH-like instability is responsible for the frequent outbursts. 


5.4. Angular momentum transport by waves 

In Fig. 1^ we visualized the total AM flux in the disk, which con¬ 
sists of the advective part, the viscous part in the disk, and the 
Reynolds part (see Eq. ([T7|l). It is, however, of great interest to 
see how much the instabilities contribute to the total AM flux 
and to deduce the efficiency of this transport. As we have de¬ 
scribed in Sect. the instability in the BL leads to the devel¬ 
opment of acoustic waves that are launched from the interface 
and propagate into the disk and the star. As a consequence, AM 
transport in our models is dominated by acoustic radiation, i.e. 
sound waves carrying AM. We do not observe any turbulence in 
our simulations and hence rule out turbulent stresses as a ma¬ 
jor contribution to the AM flux. Unlike the turbulent stresses, 
the AM transport by waves is a non-local process. Waves can 
carry AM of either sign from one point in the domain, e.g. the 
interface in the BL where they are launched to another point in 
the disk where they deposit the AM into the fluid through wave 
dissipation. Thus, the long-range exchange of angular momen¬ 


tum between regions that are otherwise uncorrelated is possible 
through the wave transport. This is a major difference from the 
cr-models which yield a local effective viscosity due to a local 
AM exchange phenomenon, such as turbulence. Therefore, it is 
questionable whether the processes in the BL can be described 
by a local prescription like the a-model. 

We now analyze the efficiency of the wave mediated AM 
transport and consider Fig. [^for this purpose. It shows the di¬ 
mensionless Reynolds stresses q-re, which has been calculated 
according to Eq. ( [I4) i, as a function of radius and time. We can 
apply this description to the AM transport since waves trans¬ 
port angular momentum exclusively by stresses. The Reynolds 
stress is evidently very strong with peak values up to the order 
of unity and it is largest at the top of the BL, at around r = 1.2, 
which is very interesting since in that region Q. reaches its max¬ 
imum. Thus, our simulations reveal that the stresses are largest 
around the very point where they would vanish according to a 
conventional a-viscosity model. Strictly speaking, this renders a 
local prescription that depends on the shear of the angular veloc¬ 
ity inapplicable for the BL. Eigure also shows that the AM 
transport due to the waves includes a large part of the domain, 
1.05 < r < 1.4, and spreads out farther into the disk during out¬ 
bursts. This reinforces the point that the AM transport ultimately 
generated by the shear in the BL is not confined to the BL, but 
in fact ranges over long distances. As we saw earlier, the activity 
rises during outbursts and therefore ctre is enhanced during these 
periods as well. The value of aRE is almost exclusively negative, 
meaning that angular momentum is transported toward the star. 
This behavior is characteristic of the lower branch of the acoustic 
modes ( Belyaev et al.|2013a i. We refrain from presenting plots 
of the parameter ay that can be calculated according to Eq. 
because it is artificial to boil the non-local AM transport down 
to a local parameter. However, for the sake of completeness, it 
should be mentioned that Eq. ( [T5] l yields values on the order of 
10 “^ - 1 for our models. 

Einally, we investigate how the strong AM transport in and 
around the BL affects the width of the region where the az- 
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Time [orbits] 

Fig. 15. Time-radius image of ttRe spanning several hundred orbits. We note that the values of ctrc have been multiplied by (-1) in order to allow 
for the logarithmic scale. For each point in time the values have been time average over ten orbits. The vertical dashed lines denote the locations 
of the outbursts mentioned earlier in the text. 


imuthal velocity rises from stellar rotation to the Keplerian ro¬ 
tation in the disk. Figure [^depicts the BL thickness as a func¬ 
tion of time for all four simulations. The thickness of the BL is 
dehned to be the region in which 

0.02 < u^ir)luii(r) < 0.8. (26) 

The width of the BL increases very rapidly and grows in size 
by a factor of 2 at the beginning of the simulation, i.e. when the 
sonic instability has grown to establish an efficient AM trans¬ 
port. After the initial fast growth the width of the BL increases 
on average linearly over the course of approximately 500 orbits 
until it settles down at a value of ~ Q.2R,. We have run additional 
simulations in order to test the robustness of the BL thickness of 
20% of the stellar radius. In these simulations, the initial dis¬ 
tance of the velocity drop from the inner edge of the simulation 
domain has been extended by altering both the parameter 'F (see 
Sect. |3.2| l and the point where the computational domain begins. 
We thereby test whether the thickness is affected by the BL run¬ 
ning up against the inner edge of the domain. However, these 
runs are in perfect agreement with the simulations presented here 
and saturate at a BL width of ~ Q.2Rt as well. 

There are several spikes in the temporal evolution of the BL 
width for all three simulations that are clearly associated with the 
frequent outbursts. During the high activity state the BL broad¬ 
ens considerably in a very short period of time and typically re¬ 
sides at a larger width in the aftermath of the outburst. Because 
of the high resolution and the resulting long simulation times 
we can follow the highly resolved simulation only for about 200 
orbits. During that time it behaves exactly like the other simula¬ 
tions with the exception that the width of the BL increases much 
more quickly. This is in good agreement with the fact that the 
growth rate of the sonic instability depends on the resolution in 
such a way that it is faster when the the resolution is higher (see 
Sect.|^. Since the BL widens to a large extent, it extends very far 
into the disk and into the star. The latter means that a consistent 
simultaneous treatment of the star is required for numerical sim¬ 
ulations involving large simulation times. We will discuss this 
point in detail in the next section. 

5.5. The long-term evolution of our simulations 

We ran and monitored our reference simulation for over 2000 
orbits (see Fig. |^, and the instability and the periodic outbursts 
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do not die out during that time. There is, however, a change in 
the pattern of the acoustic waves that occurs around orbit 1200, 
which coincides precisely with the moment when the base of 
the BL reaches the inner boundary of the simulation domain. 
We consider this a problematic moment, since it is not possible 
to treat the star correctly in 2D r-tp simulations. The geometry 
of the coordinate system inherently gives the star a cylindrical 
shape in our simulations. In Sect. |3.2| we describe how we treat 
the inner boundary in order to approximate the beginning of the 
star and that the radial velocity is set to outflow with a certain ve¬ 
locity at this point. Therefore, we might expect unphysical and 
artihcial processes once the vortices of the acoustic modes reach 
this point because the radial velocity is fixed to a certain value 
and an inflow in the domain is not supported. Furthermore, we 
force the azimuthal velocity to zero at the inner boundary. How¬ 
ever, the acoustic modes seem to be quite robust since they do not 
die out despite the difficulties with the artihcial inner boundary. 
The change in the pattern at orbit 1200 might be due to a switch¬ 
ing of the modes to the upper branch, since the spiral arms of the 
waves are unwinding considerably and changing from leading to 
trailing. Another 1000 orbits later the process is reversed again 
and the system switches back to the lower branch. 

However, in order to capture the physics correctly at these 
late stages in time and to still make reliable statements once the 
BL extends all the way to the surface of the star or even into 
the star, it is necessary to attach a realistic model of the star to 
the disk. We consider, for instance, a wave launched in the BL 
that propagates toward the star. It can penetrate deep into the star 
before it is stopped by the steep density gradient and dissipates 
the angular momentum somewhere in the star. On the other hand, 
many problems arise when the star is treated consistently. The 
problem must then be simulated in spherical coordinates, i.e. 3D 
simulations are necessary, and the density varies considerably 
from the interior of the star to the disk. Simulations of this kind 
are very demanding from a computational point of view, but also 
very tempting, and we plan to investigate this problem in the 
future. 

6. Summary and conclusion 

We have performed highly resolved, 2D hydrodynamical simu¬ 
lations of the BL surrounding a young star in order to investigate 
the AM transport driven by instabilities. Extending previous sim- 
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ulations, we have a net mass flow through the disk, we utilized 
a realistic equation of state, and (within the framework of one- 
temperature approximation) we included full radiative transfer 
in the disk plane as well as vertical cooling through an realis¬ 
tic estimate of the vertical optical depth, thereby employing a 
quasi 3D radiation transport. Moreover, those simulations were 
started from state-of-the-art ID models of the BL and thus com¬ 
prise realistic profiles for the density and temperature. We ran 
the simulations for over two-thousand orbits in order to study 
the long-term behavior. 

Our simulations confirm that the supersonic velocity drop in 
the BL is indeed prone to the sonic instability, a kind of super¬ 
sonic shear layer instability whose subsonic counterpart is the 
Kelvin-Helmholty instability (Belyaev & Rafikov 2012| l. Shortly 
after starting the simulations, the sonic instability sets in and 
starts to grow rapidly for about 15 orbits. It then saturates and 
the BL is dominated by acoustic waves that propagate both into 
the disk and into the star. These acoustic modes manifest them¬ 
selves in three distinct patterns that can be related to the three 
branches of the instability of a plane parallel vortex sheet and do 
not die out for the whole simulation time. 

Additionally, the system repeatedly undergoes outbursts 
where the wave activity as well as the AM and mass transport 
increase considerably. We argue that these outburst are likely 
triggered by a secondary KH instability that develops in the flat 
region of the azimuthal velocity profile due to several changes of 
sign of the vorticity. Since the effective temperature also shows 
strong variations during outbursts, they may play an important 
role in explaining variations in the light curve, such as FU Or- 
outbursts or DNOs and QPOs in the case of BLs around white 
dwarfs. It is tempting to associate these phenomena with the out¬ 
bursts we observe in our simulations; however, this topic must be 
looked into further in order to draw reliable conclusions. Such an 
investigation would certainly involve the study of how the out¬ 
burst reacts to a change in parameters, among other things. It is 
also vital to perform 3D simulations on this question since there 
might be a dependence on dimensionality ( [Belyaev et al.|2013^ . 

Our main objective in this work was the clarification of the 
AM transport in unmagnetized astrophysical BLs. We found that 
the acoustic modes indeed transport AM through the BL and the 
disk, and that this mechanism is even highly efficient, reaching 
values for the dimensionless Reynolds stress aRe of ~ 0.01 in 
quiet states and up to unity in high activity states (see Fig. 
Along with AM, mass is also transported efficiently through the 
domain (see Fig.[7]l. The BL reacts to the enhanced transport by 
widening considerably until it reaches a thickness of about Q.2Rt 
(see Fig.fT^. 

One of the most important elements of wave mediated AM 
transport is the intrinsic non-locality of this process, i.e. waves 
can extract AM from one point in the disk and release it in an¬ 
other farther off. There are several implications that emerge from 
this: Thus far, or-models for the parametrization of the viscosity 
have been used both in the disk and in the BL, possibly with 
some corrections for the BL. This prescription is justified in the 
disk where, under certain conditions, the system is susceptible 
to the MRI. The MRI leads to turbulence in the gas and angular 
momentum is transported via turbulent stresses. This is a local 
process that also depends on the shearing at this point and a lo¬ 
cal effective viscosity can be derived by the model of jShakura] 
& Sunyaev (|1973|). Furthermore, it allows for efficient mixing of 


the gas through the turbulence. In contrast, the AM transport we 
have investigated in this work does not depend on the local shear¬ 
ing, for instance, but is a non-local mechanism, instead. There¬ 
fore, it is doubtful whether a conventional a-model is applicable 


in the BL, and it will be hard to give a simple parametrization for 
the AM transport in the BL at all. However, aRe features a rather 
smooth behavior and might be utilized for a simplified picture of 
the stresses in the BL. In addition, the mixing in the BL might 
not be as pronounced as in the disk, which could have important 
implications concerning the spectrum of the BL. 

Connected with the notion of non-local AM transport is the 
fact that the kinetic energy, which resides in the gas shortly be¬ 
fore it is decelerated to stellar rotation, does not need to be re¬ 
leased as locally as was previously assumed. Waves also trans¬ 
port energy and release it to the fluid where they are damped or 
dissipate. This might have important implications for the struc¬ 
ture and observational appearance of the BL. We consider, for 


instance, the alternative BL theory of the spreading layer (Inog- 
amov & Sunyaev||l999| 2010|l where the gas is not decelerated 


in the disk midplane but rather in two belts in the northern and 
southern hemisphere of the star. It is difficult to bring a local vis¬ 
cosity model into accordance with this theory, since the friction 
on the surface of the neutron star is far too small. The non-local 
AM transport described in this work could, however, shed new 
light upon the theory of the spreading layer. 

Simulations that investigate both the vertical structure and 
the non-axisymmetric properties can only be performed in 3D, 
however. This is a vital extension that should be taken as the next 
step in order to bring the models closer to reality. In a 3D spheri¬ 
cal geometry, it would also be possible to simultaneously model 
the star, a point that has been illustrated in Sect. |5.5| Another ma¬ 
jor extension is the inclusion of magnetic fields, i.e. performing 
magnetohydrodynamical simulations of the BL. It would then be 
possible to relinquish any viscosity prescription whatsoever and 
directly simulate the MRI in the disk and the acoustic modes in 
the BL. There might be important interactions between the tur¬ 
bulent disc and the BL modes and new wave branches might add 
to the three discussed here. However, the extensions mentioned 
above will greatly increase the computation time and possibly 
require more efficient codes and faster hardware. 
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